arXiv: 1502.00784vl [physics.ao-ph] 3 Feb 2015 


The representation of bonndary currents in a finite element 

shallow water model 

Peter D. Diiben*and Peter Korn^ 

Max Planck Institute for Meteorology, 

Bundesstrasse 55, 20146 Hamburg, Germany 

February 4, 2015 


Abstract 

We evaluate the influence of local resolution, eddy viscosity, coastline structure, and boundary 
conditions on the numerical representation of boundary currents in a finite element shallow-water 
model. The use of finite element discretization methods offers a higher flexibility compared to finite 
difference and finite volume methods, that are mainly used in previous publications. This is true for 
the geometry of the coast lines and for the realization of boundary conditions. For our investigations 
we simulate steady separation of western boundary currents from idealized and realistic coast lines. 
The use of grid refinement allows a detailed investigation of boundary separation at reasonable 
numerical cost. 


1 Introduction 

The mechanisms that lead to separation of boundary currents in the ocean are poorly understood. Nu¬ 
merical models can provide only a small contribution to a better understanding of boundary separation, 
since the representation of the coast line in today’s ocean models differs in essential properties from the 
coast line in real-world oceans, mainly due to the coarse resolution, the high viscosity, and the imprint 
of the used grid structure. This paper investigates boundary separation in finite element models. It is 
possible that finite element discretization methods will be a common choice to build up future ocean 
models (see for example the model development in |DKS041 |PPG~*~0^ ICLDLO^ IDKA12] 1. 

For the separation of boundary currents in the ocean, such as the Gulf Stream, there are many 
possible mechanisms that might influence the position of the separation point such as a change of the 
direction in the wind field, a potential vorticity crisis, an adverse pressure gradient, boundary conditions, 
a collision with another western boundary current, interactions with the deep western boundary current, 
the coast line geometry, the bottom topography or eddy-topography interactions (see for example |Sto48[ 
IGIY871 ICCMl IEM9211HMG921 l()GP971 ITMOOl IMMOhl IGMn8| and the references therein). 

While the Gulf Stream tends to overshoot the separation point of the real world in standard numerical 
ocean models, state-of-the-art high-resolution model simulations, with a grid resolution of 1/10° or 
higher, obtain an improved representation of Gulf Stream separation |BHS071 IGGOlj . However, high- 
resolution does not guarantee a proper representation of the Gulf Stream, and the separation point 
remains sensitive to changes in the model setup, such as changes in viscosity parameterization [BHS07] . 
The choice of boundary conditions is also known to have a significant influence on the separation behavior 
|Den931lHMGM] . 

The two discretization methods that are widely used in today’s ocean models are the finite difference 
and the finite volume method. The finite difference method offers only a poor representation of the coast 
line. To introduce a coast line into a finite difference model, grid points on land are typically removed 
from a fixed grid. Structured longitude/latitude grids which are typically used allow an angle of 0 or 
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90 degrees between neighboring grid edges along the boundary. This leads to staircase pattern along 
the coast line. Furthermore, due to the staggering of the velocity components, the effective boundary 
conditions can be dependent on the angle between the coast line and the coordinate axis of the numerical 
grid (see |AM98| for the analysis on an Arakawa C-grid and B-grid). Finite volume and finite element 
methods offer higher geometric flexibility. In finite volume methods, the velocity field is represented 
as one dimensional vector perpendicular to grid edges. In finite element methods, the velocity field is 
typically defined as two-dimensional vector quantity all along the coast line. Both methods allow the 
use of boundary conforming grid generators, in which the vertices at the boundary of a grid are aligned 
to the coast line. Despite the improved coast line representation, a detailed analysis of the properties of 
boundary currents and boundary separation has not been done for finite element models with realistic 
coast lines for ocean models. 

We study the numerical representation of western boundary currents in a finite element model and 
compare the results to finite difference simulations from the literature. The used finite element model uses 
a discontinuous linear representation for velocity and a continuous second order representation for height. 
It was developed for the particular use in atmosphere and ocean model and fulfills the Ladyzhenskaya- 
Babuska-Brezzi-condition - which is a necessary condition for convergence in finite element modeling 
- and is able to represent the geostrophic balance at the same time [CHP091ICHPR09] . We simulate 
the separation of steady western boundary currents from idealized coast lines, and coast lines as used 
in ocean models. We vary the resolution, the eddy viscosity, the grid structure, the coast line, the 
alignment between the velocity components and the coast line, and between no-slip and free-slip boundary 
conditions. We evaluate the influence of these properties on boundary currents, and boundary separation. 
The test cases studied in this publication do not fundamentally differ from test setups used in publications 
such as [Den93] . [HMG92) . or [OCP97) for simulations with finite difference models with vorticity as 
prognostic quantity. The main differences are that we use a finite element model, and velocity and 
height as prognostic quantities. 

In section two, we give a very short description of the model setup, including the shallow-water 
equations, the discretization in space and time, and grid refinement. In section three, we introduce the 
test cases and present the numerical results. In section four, we discuss the results. 

2 Model setup 

This section provides a brief introduction to the functionality of the used model, including the shallow- 
water equations, the discretization in space and time, and the used grids. A detailed description of the 
model setup can be found in [DKAI2) . 

2.1 The viscous shallow-water equations 

Our finite element model simulates the viscous shallow-water equations in non-conservative form 

I r® 

dtu + u • Vu + fkx u + gVh —— V ■ (Hv'Vu) = — — y/u, (1) 

H H 

dth + y ■ (Hu) = 0, 

where u is the two dimensional velocity vector, / is the Coriolis parameter, k is the vertical unit vector, 
g is the gravitational acceleration, v is the eddy viscosity, r'* is the surface wind forcing, 7 / is the bottom 
friction coefficient, h is the surface elevation and H is the height of the fluid column given hy H = h — ht-, 
where ht is the bathymetry. The prognostic variables are the surface elevation and the velocity. 

The used model can run with either free-slip (u • n = 0, and (?nU = 0 on the boundary 9D), or no-slip 
boundary conditions (u = 0 on dQ). We apply a weak representation of no-slip boundary conditions of 
zero tangential velocity. To this end, we add the penalty term —cru to right hand side of equation [T] for 
all velocities along the boundary that pushes the tangential velocity along the boundary towards zero. 
O’ is a constant that needs to be adjusted experimentally. All other boundary conditions are realized as 
strong boundary condition by adjusting the corresponding numerical fluxes through the boundary. 
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2.2 Discretization in space and time 

Following the typical finite element approach, we expand the physical fields into sets of basis functions 
Ni and Mi 


Nu Nh 

u = viiNi and h = hiMi. 

i=l i=l 

We use a PP^P 2 finite element approach to discretize the equations. This means that we employ 
discontinuous linear Lagrange polynomials for the representation of the velocity field (Ni), and globally 
continuous quadratic Lagrange polynomials for the representation of the height field {Mi). Each trian¬ 
gular cell has three degrees of freedom for each component of velocity located at the vertices of the cells, 
and six degrees of freedom for the height field located at the vertices and edges. While the degrees of 
freedom of the height field are shared with the surrounding cells, the degrees of freedom of the velocity 
field belong to a specific cell, which leads to a discontinuous representation. 

Time integration is performed by an explicit three-level Adams-Bashforth method. The equation 

dt-ip = R{ip), 

where R denotes the right-hand-side of the system, and '0 is the vector of prognostic variables, is 
discretized in time by 


0-+1 = 0- -f At (^§i?(0") - ^i?(0"-') + ^^(^"■')) . 
where 0* is the vector of state variables at the i-th time step. 

2.3 Grids and grid refinement 

We use two types of standard grids on which refinement is performed. The first type of grids are 
structured triangular grids that provide a uniform coverage of the longitude/latitude space. The grids 
are derived from rectangular grids by bisecting each rectangular into two triangles. The second type of 
grids are icosahedral geodesic grids that provide a quasi-uniform coverage of the sphere [BF85] . We use 
static h-refinement to refine the interesting area around the coast line. In h-refinement, new grid points 
are introduced to the grid, to increase the model resolution in regions of specific interest. The influence 
of grid refinement to the model solution is investigated in |DK14j . 

3 Numerical tests and results 

In this section, we present the numerical results for two test cases. We will start with an idealized case 
that uses straight boundaries, simulating a steady wind driven ocean gyre with a western boundary 
current that separates at the corner of an obstacle. In the second test, we study a more realistic setup to 
investigate coast lines as used in ocean models simulating a steady wind driven circulation in an Atlantic 
shaped basin. 

3.1 Ocean gyre with idealized coast lines 

We study an ocean gyre setup on the northern hemisphere. The gyre is driven by a wind forcing in 
clockwise direction. Due to the change of the Coriolis parameter in the meridional direction, the gyre is 
intensified towards the western boundary, and a western boundary current develops |Sto481 IPed96] . The 
current is separating from the coast at the edge of a rectangular obstacle. The setup is chosen to be as 
close as possible to the setup used in [Den98) . Here, Dengg investigated boundary separation in a model 
that simulates the barotropic vorticity equation. 

We perform model runs on a triangular grid, which is structured in longitude/latitude space. We use 
refinement to increase the resolution in the vicinity of the boundary (see Figure [T]). A grid edge has a 
length of about 1.6° in the coarsest and 0.4° in the finest part of the grid, this corresponds to about 172 
and 43 km at the southern boundary of the domain. 
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Figure 1: Vertices of the grid with two refinement levels used for the idealized coast line test. The red 
line marks the coast line. 


While the meridional wind forcing is zero, the zonal wind forcing is set to 

'7r(6»- 15°)' 


"Ta = T-o 


10 ' 


cos 


40° 


where 9 is the latitude. The bottom friction coefficient 7 / is set to 10 ° s and tq is 0.28 m^s ^ 
The height field is initialized with a constant water depth of 1000 m; the initial velocity is set to zero. 


a, no—slip 1 lev. ref. b, no—slip 2 lev. ref. c, free—slip 1 lev. ref. 



999.85 1000 1000.15 1000.3 1000.45 


[m] 


Figure 2: Equilibrium height field of the idealized coast line test. While the runs a, 6 , c, and d were 
simulated with v = 3000.0 runs e and/ were simulated with ly = 10000.0 The simulations 

are performed on the grid plotted in Figure [1] either without refinement (e), with one refinement level 
(a and c), or with two refinement levels (&, d, and /). 

Figured] shows the equilibrated height field. For all tests, the Munk layer at the western boundary 
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is represented smoothly. 

The simulations in [Den93) show no boundary separation when free-slip boundary conditions are used. 
Nevertheless, in the present simulations the boundary flows separate for free-slip and no-slip boundary 
conditions. The equilibrated fields have a different shape for the two different boundary conditions 
(compare b and d). While resolution does not play an important role for boundary separation (compare 
a with b, c with d, and e with /), changes in viscosity have a strong impact (compare d and /). 
The simulated flows have Reynoldsnumbers between 10 and 100 along the boundary. The results are 
qualitatively the same when simulations are performed with stronger wind forcings, and therefore with 
higher Reynoldsnumbers. We do not show these results since the resulting equilibrium fields are unsteady 
and it is much more difficult to compare them. 



Longitude 


Longitude 


Longitude 


Longitude 
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1000 


1000.1 1000.2 1000.3 1000.4 1000.5 
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Figure 3: Equilibrium height field of the idealized coast line model runs with no-slip (top row) and 
free-slip (bottom row) boundaries conditions on different grids. Results should be compared with the 
model runs a and c in Figure [5] 

The simulations in Figure [3] are performed with the same setup as for the previous idealized coast 
line tests a and c in Figure[2]with i/ = 3000 tq = 0.28m^s~^. However, different grids have been 

used. Simulations a and d were performed on a one level refined icosahedral grid (see Figure S]). The 
coast line is represented by an unstructured pattern. The simulations b and f were performed on a grid 
with a changed longitude of the model domain compared to the standard grid. A change of the longitude 
is changing the alignment of the two components of the velocity field u and v with the coast line. In the 
model runs c, d, g and h the arrangement of the triangles in the structured grid was changed (Figure [5]). 
The change of the grid structure leads to a zig-zag pattern of the coast line in model run d and h. 

While all simulations with no-slip boundary conditions result in a fairly similar gyre structure, this 
is different for the free-slip simulations. The simulations with free-slip boundary conditions and unstruc¬ 
tured or zig-zag meridional coast line (e and h) appear to be similar to the no-slip simulations. A likely 
reason for this similarity is a shift of the boundary flow into the interior of the domain that is caused 
by the boundary conditions in the no-slip case and by the abrupt changes of the direction of the coast 
line in simulations e and h. A similar behavior can be found for simulations with ocean models based 
on an Arakawa C-grid finite difference scheme that change the boundary conditions from free-slip to 
no-slip when zig-zag coast lines are considered (see [AMQSj b but the mechanism is different. In the finite 
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Figure 4: Vertices of the refined grid built from an icosahedral grid, used in model run a in Figure [31 



Figure 5: Structure of the triangular grids used for the standard runs (left), simulations c and g (middle), 
and simulations d and h in Figure |3] (right), with indicated coast line at the western boundary. 


element model, the representation of the boundary conditions stays the same when the alignment of the 
velocity components with the boundary change. This is confirmed by the simulations b and f that show 
almost identical flow pattern compared to simulations a and c in Figure O 

While the separation behavior in the free-slip simulation g in Figure [3] shows differences to the 
reference simulation a in Figure [21 the separation behavior of the no-slip simulation c is very similar to 
the reference simulation. The slight change of the coast line is able to change the solution of the free-slip 
simulation, while this is not true for the no-slip simulation. 


3.2 Irregular coast lines - An Atlantic shaped ocean domain 

In this section, we investigate a more realistic application to model the ocean. We simulate an ocean 
basin which is shaped like the Atlantic ocean and offers a realistic representation of the coast line. The 
domain is cut at the equator and at 58° North. We simulate on real-world topography, but water depth is 
cut at 1000 m. An artificial wind forcing that is balanced by bottom friction induces a steady circulation. 
The used numerical grid is plotted in Figure [51 The grid is refined at the western boundary and has a 
typical edge length of 120 km in the coarse, and 60 km in the fine part of the grid. In the refined area 
along the coast line there are always two neighbored grid edges that are aligned with each other. 

Simulations are initialized with zero surface elevation and zero velocity. The zonal wind forcing is 
given by 


-To ■ 10-3 ■ cos (4 • 0) 

0 


if 0 < 45° 

if 0> 45°, 


the meridional wind forcing is zero. The bottom friction coefficient 7 / is set to 10 ® s and tq is 
3.0 
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Figure 6: Vertices of the refined grid used for the Atlantic test. The red line marks the coast line. 



994 995.5 997 998.5 1000 1001.5 


[m] 

Figure 7: Equilibrium height field in the Atlantic model runs, after 140 days. While run a and b were 
simulated with i/ = 6655.0 run c and d were simulated with v = 53240.0 
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Figure [7] shows the equilibrium height field of the performed model runs. 

The model runs a, b and c are performed on the grid plotted in Figure IHl Model run d is using 
the same grid without refinement of the western boundary. As for the simulation with an unstructured 
representation of the idealized coast line (a in Figure [3]) the unstructured realistic coast line seems to 
make the used boundary condition play only a minor role - the height field along the western boundary 
differs much more for different values of eddy viscosity (compare a with c) than for no-slip and free-slip 
boundary conditions (compare a and b). A change in resolution leads to minor changes (compare c and 
d), although it limits the smallest possible value for eddy viscosity (see also |DK14] 1. 

4 Discussion of the results 

Although finite element methods provide an improved coast line representation compared to finite dif¬ 
ference methods, our investigations show that the representation of the coast line and the boundary 
conditions is still not satisfying. Changes of the grid structure can lead to changes of the separation 
behavior (see Figure [S]). 

Our tests on the influence of resolution and eddy viscosity show that steady western boundary cur¬ 
rents are not strongly affected by changes in resolution, as long as the Munk layer is resolved properly. 
However, a higher resolution allows the use of a smaller eddy viscosity, which can change the model 
results significantly. To this end, grid refinement can be used to increase the local resolution in the 
Munk layer. 

The model results change strongly between free-slip and no-slip boundary conditions, when idealized 
straight coast lines are simulated (subsection [OJ- This is heavily discussed in the literature (see |Den93) 
as one example). On the other hand, the model results change only slightly with the boundary conditions 
for coast lines as used in ocean models (subsection 13.2|) . In simulations with free-slip boundaries and 
zig-zag or unstructured coast lines, the flow is shifted towards the interior of the domain, due to the 
rapid changes of the direction of the coast line. The results look similar to no-slip model runs, in which 
the flow is shifted into the interior of the domain via the boundary conditions Isubsection 13.11) . 

In contrast to finite difference models with the vorticity as prognostic quantity |Den93] . we obtain 
separation for free-slip boundary conditions, using a finite element model with velocity and height as 
prognostic quantities. We do obtain premature separation for no-slip, but not for free-slip boundary 
conditions Isubsection 13.111 . This result is consistent with results of finite difference models. 

Although finite element methods offer an improved coast line representation compared to finite dif¬ 
ference methods, the representation of boundary flows remains dependent on the pattern of the coast 
line, which is - for today’s ocean models - very much dependent on the resolution, and not satisfying. 
Small changes of the grid structure can lead to changes in the separation behavior. 
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